---
title: "Camelid size"
author: "Juliana Rubinatto Serrano"
date: "2024-01-29"
output: html_document
---

Este compendium fue utilizado por el análisis del trabajo por:

deFrance, S. D y Juliana Rubinatto Serrano (2025) "Camélidos Tihuanacos en la Colonia de Moquegua al Sur del Perú: Variaciones de Tamaño y Morfotipos". *Chungara: Humanos y camélidos: Interacciones sociales e historia evolutiva.*

# Datos y "*packages*" utilizados

```{r setup}
library(tidyr)
library(dplyr)
library(reshape2)
library(zoolog)
library(ggplot2)
library(ggpubr)
library(rstatix)
```

```{r data wrangling}
data <- read.csv("data/data.csv")

```

# Análisis exploratorio


```{r EDA graphics}

data$species <- factor(data$species,levels = c("Moquegua camélido (n = 298)", "Le Neün et al. (2023) camélido arqueologico (n = 125)", "Alpaca moderna (n = 22)", "Guanaco moderno (n = 179)", "Llama moderna (n = 109)", "Vicuña moderna (n = 109)"))

# Por todas las muestras (usado para generar Figura 3 en deFrance and Rubinatto Serrano (2025))
plot_by_sample <- ggplot(data, aes(x=H, y=Bp, 
                                 shape = species, 
                                 color = species)) + 
  geom_point(size = 10) + 
  annotate(geom = "text",x=13, y=10, label="Pequeño", size = 15) +
  annotate(geom = "text",x=25, y=10, label="Grande", size = 15) +
  scale_color_manual(values=c("#7851A9", "#FF7900", "#DF30FF",
                              "#48ADDE", "#A5DE00", "#FF9111"),
                     labels = c("Moquegua camélido (n = 298)", 
                                "Le Neün et al. (2023) \ncamélido arqueologico (n = 125)",
                                "Alpaca moderna (n = 22)", 
                                "Guanaco moderno (n = 179)", 
                                "Llama moderna (n = 109)", 
                                "Vicuña moderna (n = 109)")) +
  scale_shape_manual(values=c(19, 20, 18, 18, 18, 18),
                     labels = c("Moquegua camélido (n = 298)", 
                                "Le Neün et al. (2023) \ncamélido arqueologico (n = 125)",
                                "Alpaca moderna (n = 22)", 
                                "Guanaco moderno (n = 179)", 
                                "Llama moderna (n = 109)", 
                                "Vicuña moderna (n = 109)"))
  
plot_by_sample + labs(x = "V3 (mm)", y = "V2 (mm)", color = "", shape = "") + 
  theme(legend.position = "bottom",
        axis.title = element_text(size = 35), 
        plot.title = element_text(size = 35), 
        legend.text = element_text(size = 40))+
  geom_abline(intercept = 45.5, slope = -1.5, color="red", 
                 linetype="dotdash", size=.5) +
  geom_abline(intercept = 42, slope = -1.5, color="red", 
                 linetype="twodash", size=.5)


# Por los sitios Tihuanaco en Moquegua (usado para generar Figura 2 en deFrance and Rubinatto Serrano (2025))

Moquegua <- data %>% filter(group == "Moquegua (n = 298)")
long_data_Moquegua <- Moquegua %>% select(Sitio, Bp, H) %>% 
  gather(key = "variable", value = "value", -Sitio) #long format
long_data_Moquegua_GL <- Moquegua %>% select(Sitio, GL) %>% 
  gather(key = "variable", value = "value", -Sitio)

species <- data %>% filter(group != "Le Neün et al. (2023) arqueologica (n = 125)")

species$species <- factor(species$species,levels = c("camélido", "alpaca", 
                                                     "vicuña", "llama", "guanaco"))
long_data_species <- species %>% select(species, Bp, H) %>% 
  gather(key = "variable", value = "value", -species) #long format
long_data_speciesGL <- species %>% select(species, GL) %>% 
  gather(key = "variable", value = "value", -species) #long format


plot_by_site <- ggplot(Moquegua, aes(x=H, y=Bp, color= Sitio)) + 
  geom_point(size = 10) + 
  annotate(geom = "text",x=13, y=10, label="Pequeño", size = 15) +
  annotate(geom = "text",x=25, y=10, label="Grande", size = 15) +
  scale_color_manual(values=c("#DF30FF","#48ADDE", "#A5DE00", "#FF9111"),
                     labels=c("M1 - Chen Chen\n (n = 22)", "M10 - Omo\n (n = 71)",
                              "M12 - Omo\n (n = 8)","M43 - Rio Muerto\n (n = 197)"))

plot_by_site + labs(x = "V3 (mm)", y = "V2 (mm)") +
  theme(legend.position = "bottom", axis.title = element_text(size = 35),
        legend.title = element_blank(),
        plot.title = element_text(size = 35), legend.text = element_text(size = 40)) +
  geom_abline(intercept = 45.5, slope = -1.5, color="red", 
                 linetype="dotdash", size=.5) +
  geom_abline(intercept = 42, slope = -1.5, color="red", 
                 linetype="twodash", size=.5)


# Diagrama de caja por los sitios y todas las muestras (usado para generar Figura 4 en deFrance and Rubinatto Serrano (2025))

##Sitios
boxplot1 <- ggplot(long_data_Moquegua, aes(x = variable, y=value, color = Sitio)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#DF30FF","#48ADDE", "#A5DE00", "#FF9111"),
                     labels=c("M1 - Chen Chen\n (n = 22)", "M10 - Omo\n (n = 71)",
                              "M12 - Omo\n (n = 8)","M43 - Rio Muerto\n (n = 197)")) + 
  labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
  theme(legend.position = "none", legend.title = element_blank(),
        axis.title = element_text(size = 35),
        plot.title = element_text(size = 35), legend.text = element_text(size = 40))

boxplot2 <- ggplot(long_data_Moquegua_GL, aes(x = variable, y=value, color = Sitio)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#DF30FF","#48ADDE", "#A5DE00", "#FF9111"),
                     labels=c("M1 - Chen Chen\n (n = 22)", "M10 - Omo\n (n = 71)",
                              "M12 - Omo\n (n = 8)","M43 - Rio Muerto\n (n = 197)")) + 
  labs(x = "GL = V1", y = "") +
  theme(legend.position = "none", axis.title = element_text(size = 35))

##Muestras
long_data <- data %>% select(group, Bp, H) %>% 
  gather(key = "variable", value = "value", -group) #long format
long_dataGL <- data %>% select(group, GL) %>% 
  gather(key = "variable", value = "value", -group)

boxplot3 <- ggplot(long_data, aes(x = variable, y=value, color = group)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#43984A", "#FF7900", "#7851A9"),
                      labels=c("Le Neün et al. (2023)\narqueologica (n = 125)", 
                              "Le Neün et al. (2023)\nmoderna (n = 419)",
                              "Moquegua \n(n = 298)")) + 
  labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
  theme(legend.position = "none", legend.title = element_blank(),
        axis.title = element_text(size = 35), 
        plot.title = element_text(size = 35), legend.text = element_text(size = 40))

boxplot4 <- ggplot(long_dataGL, aes(x = variable, y=value, color = group)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#43984A", "#FF7900", "#7851A9"),
                      labels=c("Le Neün et al. (2023)\narqueologica (n = 125)", 
                              "Le Neün et al. (2023)\nmoderna (n = 419)",
                              "Moquegua \n(n = 298)")) + 
  labs(x = "GL = V1", y = "") +
  theme(legend.position = "none", axis.title = element_text(size = 35))



###Combinar los diagramas
big1 <- ggarrange(boxplot1, boxplot2, common.legend = TRUE, legend = "bottom")
big2 <- ggarrange(boxplot3, boxplot4, common.legend = TRUE, legend = "bottom")

ggarrange(big1, big2, nrow = 2, legend = "bottom", labels ="AUTO", font.label = list(size= 40))

```

# ANOVA, MANOVA y t-tests

```{r statistical tests}
#Entre muestras
manova_samples <- manova(cbind(GL,Bp, H) ~ group, data = data) 
summary(manova_samples)
manova_samples1 <- manova(cbind(Bp, H) ~ group, data = data) 
summary(manova_samples1)

  ##media geometrica
  data_noNA <- subset(data, !is.na(mean_geo))
  anova_geo <- oneway.test(mean_geo ~ group, data = data_noNA, var.equal = TRUE)
  anova_geo
  
  GL <- data %>% select(group, GL)
  GL_noNA <- subset(data, !is.na(GL))
  res.aov <- oneway.test(GL ~ group, data = GL_noNA, var.equal = TRUE)
res.aov

  ##two-sided t-tests entre

    ###Moquegua y le Neun et al (2023) arqueologica
    archaeological <- data %>% filter(group == "Moquegua (n = 298)" | 
                                        group == "Le Neün et al. (2023) arqueologica (n = 125)")
    arqueo <- t.test(cbind(Bp, H) ~ factor(archaeological$group), data = archaeological, 
                     alternative = "less", var.equal = TRUE)
    
    ###Moquegua y le Neun et al (2023) moderna
    mixed <- data %>% filter(group == "Moquegua (n = 298)" | 
                                        group == "Le Neün et al. (2023) moderna (n = 419)")
    mixed_test <- t.test(cbind(Bp, H) ~ factor(mixed$group), data = mixed, 
                     alternative = "greater", var.equal = TRUE)
     
    ###le Neun et al (2023) moderna y arqueologica
    leNeunetal <- data %>% filter(group == "Le Neün et al. (2023) arqueologica (n = 125)" | 
                                        group == "Le Neün et al. (2023) moderna (n = 419)")
    leNeunetal_test <- t.test(cbind(Bp, H) ~ factor(leNeunetal$group), data = leNeunetal, 
                     alternative = "less", var.equal = TRUE)
    
    print(arqueo)
    print(mixed_test)
    print(leNeunetal_test)


#Sitios Tiwanaku  en Moquegua
manova_sites <- manova(cbind(GL, Bp, H) ~ Sitio, data = Moquegua) 
summary(manova_sites)
    
manova_sites1 <- manova(cbind(Bp, H) ~ Sitio, data = Moquegua) 
summary(manova_sites1)


#Delanteras y traseras
  ##delanteras
  just_fore <- data %>% filter(fore_hind == "fore")
  manova_fore <- manova(cbind(GL, Bp, H) ~ group, data = just_fore) 
  summary(manova_fore) 
  
  ##traseras
  just_hind <- data %>% filter(fore_hind == "hind")
  manova_hind <- manova(cbind(GL, Bp, H) ~ group, data = just_hind) 
  summary(manova_hind) 


#Entre muestras de Peru y Bolivia
  Bolivia <- data %>% 
    filter(Origen == "Khonkho Wankane " | Origen == "Chiripa" | Origen == "Bolivia") %>%
    mutate(country = "Bolivia")
  Peru <- data %>% filter(source != "Moquegua" & Origen == "Peru") %>% mutate(country = "Peru")
  Moquegua2 <- Moquegua %>% mutate(country = "Tiwanaku")
  
  Peru_Bolivia <- bind_rows(Bolivia, Peru, Moquegua2)
  
  manova_origin <- manova(cbind(GL, Bp, H) ~ country, data = Peru_Bolivia) 
  summary(manova_origin) 

    ##Peru 
    Peru_t <- t.test(cbind(Bp, H) ~ country, 
                     data = (filter(Peru_Bolivia, country == "Peru" | country == "Tiwanaku")), 
                     alternative = "less", var.equal = FALSE)
    print(Peru_t)
    
    ##Moquegua and Bolivia
    Bolivia_Tiwanaku <- t.test(cbind(GL, Bp, H) ~ country, 
                               data = (filter(Peru_Bolivia, 
                                              country == "Bolivia" | country == "Tiwanaku")), 
                               alternative = "greater", var.equal = FALSE)
    print(Bolivia_Tiwanaku)
    
    ##Bolivia y Peru
    Bolivia_Peru <- t.test(cbind(Bp, H) ~ country, 
                           data = (filter(Peru_Bolivia, country == "Bolivia" | country == "Peru")), 
                           alternative = "less", var.equal = FALSE)
    print(Bolivia_Peru)
```

# Visualizaciones entre muestras de Bolivia y Perú

```{r visualization regional origin}

#Diagrama de caja de muestras de Peru, Bolivia y Moquegua (usado para generar Figura 5 en deFrance and Rubinatto Serrano (2025))

Peru_Bolivia$country <- factor(Peru_Bolivia$country, 
                               levels = c("Peru", "Bolivia", "Tiwanaku"))
  Peru_Bolivia_long <- Peru_Bolivia %>% 
    select(Bp, H, country) %>% 
    gather(key = "variable", value = "value", -country)
  Peru_Bolivia_long_GL <- Peru_Bolivia %>% 
    select(GL, country) %>% 
    gather(key = "variable", value = "value", -country)
  
  boxplot5 <- ggplot(Peru_Bolivia_long, aes(x = variable, y=value, color = country)) + 
    geom_boxplot(linewidth = 3) +
    scale_color_manual(values=c("#FDB147", "#11574A" , "#7851A9"), 
                       labels=c("Le Neün et al. (2023)\nPeru (n = 136)", 
                                "Le Neün et al. (2023)\nBolivia (n = 46)",
                                "Moquegua (n = 298)")) + 
    labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
    theme(legend.position = "bottom", legend.title = element_blank(),
          axis.title = element_text(size = 35), 
          plot.title = element_text(size = 35), legend.text = element_text(size = 40))
  
  boxplot6 <- ggplot(Peru_Bolivia_long_GL, aes(x = variable, y=value, color = country)) + 
  geom_boxplot(linewidth = 3) +
   scale_color_manual(values=c("#FDB147", "#11574A" , "#7851A9"), 
                       labels=c("Le Neün et al. (2023)\nPeru (n = 136)", 
                                "Le Neün et al. (2023)\nBolivia (n = 46)",
                                "Moquegua (n = 298)")) + 
    labs(x = "GL = V1", y = "") +
    theme(legend.position = "none", axis.title = element_text(size = 35))
  
ggarrange(boxplot5, boxplot6, common.legend = TRUE, legend = "bottom")
  
#Diagrama de caja de muestras de Bolivia y Moquegua (usado para generar Figura 6 en deFrance and Rubinatto Serrano (2025))

  Bolivia_Tiwanaku_long <- Peru_Bolivia %>% 
    filter(country == "Tiwanaku" | country == "Bolivia") %>% 
    select(Bp, H, comparison) %>% 
    gather(key = "variable", value = "value", -comparison)
  Bolivia_Tiwanaku_long_GL <- Peru_Bolivia %>% 
    filter(country == "Tiwanaku" | country == "Bolivia") %>% 
    select(GL, comparison) %>% 
    gather(key = "variable", value = "value", -comparison)
  
  boxplot7 <- ggplot(Bolivia_Tiwanaku_long, aes(x = variable, y=value, color = comparison)) + 
    geom_boxplot(linewidth = 3) +
    scale_color_manual(values=c("#FDB147", "#ADB147", "#71874A" , "#7851A9")) +
    labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
    theme(legend.position = "bottom", legend.title = element_blank(),
          axis.title = element_text(size = 35), 
          plot.title = element_text(size = 35), legend.text = element_text(size = 40))
  
  boxplot8 <- ggplot(Bolivia_Tiwanaku_long_GL, aes(x = variable, y=value, color = comparison)) + 
  geom_boxplot(linewidth = 3) +
   scale_color_manual(values=c("#FDB147", "#ADB147", "#71874A" , "#7851A9")) + 
    labs(x = "GL = V1", y = "") +
    theme(legend.position = "none", axis.title = element_text(size = 35))
  
ggarrange(boxplot7, boxplot8, common.legend = TRUE, legend = "bottom")
```

# Log size index (LSI) con el zoolog "*package*"

```{r setting up for analysis}
Measurement <- data.frame(data)

#Adicionar elemento y especies por los dados 
Measurement <- select(Measurement, -14)
Measurement$EL <- "PH1" #creating the element 1st phalanx 
Measurement$Species <- "Camelid" #creating the species name Camelid

#Separar los datos por falange delantera o trasera
Moqu_fore <- filter(Measurement, group == "Moquegua (n = 298)" & fore_hind == "fore")
Moqu_hind <- filter(Measurement, group == "Moquegua (n = 298)" & fore_hind == "hind")

#THESAURUS UPDATE: Debido a que el paquete zoolog original no incluye Camelid como taxón único, se eligió arbitrariamente a Cervus para agregar Camelid al diccionario de sinónimos. Dado que estamos usando valores de referencia personalizados, tener Camelid en Cervis no afectará los resultados, pero nos permitirá usar el paquete.

taxonThesaurus <- AddToThesaurus(
  thesaurus = zoologThesaurus$taxon,
  newName = c("Camelid"), 
  category = "Cervus")
zoologThesaurus$taxon <- taxonThesaurus
WriteThesaurusSet(thesaurus = zoologThesaurus, file = "myThesaurus")

#Adicionar las mediciones del espécimen comparativo (media de las mediciones) al package
reference1_Bp <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("Bp"),
                  Standard = c("18.89"))
reference1_Bp$Standard <- as.numeric(reference1_Bp$Standard)
reference1_H <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("H"),
                  Standard = c("17.96"))
reference1_H$Standard <- as.numeric(reference1_H$Standard)
reference1_Gl <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("GL"),
                  Standard = c("61.67"))
reference1_Gl$Standard <- as.numeric(reference1_Gl$Standard)
reference1_camelid = list(reference1_Bp, reference1_H, reference1_Gl)

#Adicionar las mediciones del espécimen comparativo (falange delantera) al package
reference2_Bp <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("Bp"),
                  Standard = c("19.92"))
reference2_Bp$Standard <- as.numeric(reference2_Bp$Standard)
reference2_H <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("H"),
                  Standard = c("18.6"))
reference2_H$Standard <- as.numeric(reference2_H$Standard)
reference2_Gl <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("GL"),
                  Standard = c("65.67"))
reference2_Gl$Standard <- as.numeric(reference2_Gl$Standard)
reference2_camelid = list(reference2_Bp, reference2_H, reference2_Gl)

#Adicionar las mediciones del espécimen comparativo (falange trasera) al package
reference3_Bp <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("Bp"),
                  Standard = c("17.86"))
reference3_Bp$Standard <- as.numeric(reference3_Bp$Standard)
reference3_H <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("H"),
                  Standard = c("17.32"))
reference3_H$Standard <- as.numeric(reference3_H$Standard)
reference3_Gl <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("GL"),
                  Standard = c("57.68"))
reference3_Gl$Standard <- as.numeric(reference3_Gl$Standard)
reference3_camelid = list(reference3_Bp, reference3_H, reference3_Gl)

```

## LSI por todas las muestras

```{r reference 1, warning = FALSE}

#LSI por Bp, H y Gl todas las muestras
testLogsBp <- LogRatios(Measurement, ref = reference1_Bp)
testLogsH  <- LogRatios(Measurement, ref = reference1_H)
testLogsGl <- LogRatios(Measurement, ref = reference1_Gl)

#Combinar los resultados
Logs <- cbind(testLogsBp, testLogsGl$logGL, testLogsH$logH)
colnames(Logs)[17] <- "logGL"
colnames(Logs)[18] <- "logH"

#Formato longo
Logs_long <- Logs %>% 
    select(group, logBp, logGL, logH) %>% 
    gather(key = "variable", value = "value", -group)

#Diagrama de caja de todas las muestras (usado para generar Figura 7 en deFrance and Rubinatto Serrano (2025))

boxplot_logs <- ggplot(Logs_long, aes(x = variable, y = value, color = group)) + 
  geom_boxplot(linewidth = 2) + coord_flip() + ylim(-0.3, 0.2) +
  scale_color_manual(values=c("#43984A", "#FF7900", "#7851A9"),
                     labels=c("Le Neün et al. (2023)\narqueologica (n = 125)", 
                              "Le Neün et al. (2023)\nmoderna (n = 419)",
                              "Moquegua \n(n = 298 por Bp y H)\n(n = 251 por GL)")) 
 
boxplot_logs + labs(x = "Diferentes medidas", y = "Relación logarítmica", color = "") +
  geom_hline(yintercept = 0, colour="red", linetype = "dashed", linewidth = 1) +
  theme(legend.position = "bottom", legend.title = element_blank(),
        axis.title = element_text(size = 35), plot.title = element_text(size = 35), 
        legend.text = element_text(size = 40))
```

## LSI por las falanges delanteras y traseras

```{r reference 2 and 3}

#Falanges delanteras
testLogsBp2 <- LogRatios(Moqu_fore, ref = reference2_Bp)
testLogsH2 <- LogRatios(Moqu_fore, ref = reference2_H)
testLogsGl2 <- LogRatios(Moqu_fore, ref = reference2_Gl)
Logs_fore <- cbind(testLogsBp2, testLogsGl2$logGL, testLogsH2$logH)
colnames(Logs_fore)[17] <- "logGL"
colnames(Logs_fore)[18] <- "logH"

#Falanges traseras
testLogsBp3 <- LogRatios(Moqu_hind, ref = reference3_Bp)
testLogsH3 <- LogRatios(Moqu_hind, ref = reference3_H)
testLogsGl3 <- LogRatios(Moqu_hind, ref = reference3_Gl)
Logs_hind <- cbind(testLogsBp3, testLogsGl3$logGL, testLogsH3$logH)
colnames(Logs_hind)[17] <- "logGL"
colnames(Logs_hind)[18] <- "logH"

#Combinar los resultados
Logs_all <- bind_rows(Logs_fore, Logs_hind)
Logs_all_long <- Logs_all %>% 
  select(fore_hind, logBp, logGL, logH) %>% 
  gather(key = "variable", value = "value", -fore_hind)


# Diagrama de dispersión (usado para generar Figura 8 en deFrance and Rubinatto Serrano (2025))

supps.labs <- c("Falange delantera (n = 6)", "Falange trasera (n = 8)")
names(supps.labs) <- c("fore", "hind")

ggplot(Logs_all_long, aes(x = variable, y = value, 
                          color = fore_hind, shape = fore_hind)) +
  geom_jitter(width = 0, height = 0, size = 5) +
  facet_grid(~ fore_hind, 
             labeller = labeller(fore_hind = supps.labs)) +
  theme(legend.position="none", 
        axis.title = element_text(size = 35), 
        plot.title = element_text(size = 35)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2, color = "red", linewidth = 1) + 
  scale_color_manual(values = c("#5851A9", "#01954A"), 
                     labels = c("Falange delantera", "Falange trasera")) +
  scale_shape_manual(values = c(17, 16), 
                     labels = c("Falange delantera", "Falange trasera")) +
  labs(x = "Relaciónes logarítmicas", y = "Valores") 

```


## LSI otros huesos

```{r}
#Los datos
others <- read_excel("data/others.xlsx")
others_long <- others %>% 
  select(Element, LSI) %>%
  gather(key = "variable", value = "value", -Element)

# Diagrama de dispersión (usado para generar Figura 9 en deFrance and Rubinatto Serrano (2025))

ggplot(others_long, aes(x = variable, y = value, color = Element, shape = Element)) +
  geom_jitter(width = 0, height = 0, size = 3) +
  facet_grid(~ Element) +
  theme(legend.position="none", axis.title = element_text(size = 14),
        axis.text.x = element_text(angle = 90)) +
  coord_flip() + scale_y_continuous(breaks=seq(0, 0.08, 0.04)) +
  geom_hline(yintercept = 0, linetype = 2, color = "red", linewidth = 1) + 
  labs(x = "Relaciónes logarítmicas", y = "Values") 
```
